A general approach to quantum dynamics using a variational master equation: 
Application to phonon-damped Rabi rotations in quantum dots 
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We develop a versatile master equation approach to describe the non-equilibrium dynamics of a 
two-level system in contact with a bosonic environment, which allows for the exploration of a wide 
range of parameter regimes within a single formalism. As an experimentally relevant example, we 
apply this technique to the study of excitonic Rabi rotations in a driven quantum dot, and compare 
its predictions to the numerical Feynman integral approach. We find excellent agreement between 
the two methods across a generally difficult range of parameters. In particular, the variational master 
equation technique captures effects usually considered to be non-perturbative, such as multi-phonon 
processes and bath-induced driving renormalisation, and can give reliable results even in regimes in 
which previous master equation approaches fail. 



Introduction - The recent experimental characterisa- 
tion of exciton-phonon interactions in a coherently-driven 
semiconductor quantum dot (QD), and their interpreta- 
tion in terms of a two-level system in contact with a 
bosonic environment ^2. have demonstrated that QDs of- 
fer a natural platform in which to explore dissipative 
dynamics in the solid-state. In particular, the inter- 
play between laser-driven coherent excitonic oscillations 
and incoherent phonon-induced processes can lead to 
a rich dynamical behaviour^— For example, the phe- 
nomenon of excitation-induced dephasing has been the 
subject of intense experimental and theoretical investi- 
gation recently— _ — These studies reveal that the driven 
exciton oscillation frequency and damping rate depend 
sensitively, and non-monotonically, on the energy scales 
of both the driven QD and the bulk phonon modes. How- 
ever, the relatively large range of driving strengths, tem- 
peratures, and magnitudes of exciton-phonon couplings 
that are being explored experimentally have not so far 
been captured by a single, intuitive, quantum master 
equation, making the interpretation of results potentially 
difficult. 

Existing techniques for calculating the density matrix 
dynamics of an open quantum system include methods 
which can numerically converge to give exact results; 
one example is the quasi-adiabatic propagator path inte- 
gral (QUAPI). 12 . However, such numerical methods can 
be very computationally expensive in certain parame- 
ter regimes, and typically require a complementary ap- 
proach in order to fully appreciate the complex dynami- 
cal behaviour they predict. In contrast, quantum master 
equations naturally relate the evolution of density ma- 
trix elements with environmental parameters in an intu- 
itive way, and are often less computationally expensive. 
However, the derivation of a practical master equation 
always involves an approximation, and a new equation is 
typically derived for each distinct parameter regime. In 
the case of driven quantum dots, popular choices have 



included a weak coupling approach^^ whose validity 
relies on weak exciton-phonon coupling; and a polaron 
transform method^&i 3 - which works for stronger cou- 
pling, but which in turn introduces conditions on the 
driving strength^ 

In this work, we develop a physically-motivated, ver- 
satile and efficient master equation approach to describe 
the dynamics of a continuously-driven QD in contact 
with a phonon environment, through a combination of 
a variationally-optimized unitary transformation and the 
time-local projection operator technique. In those limits 
appropriate for either a weak-coupling or polaron treat- 
ment we find that our master equation yields similar re- 
sults, but it is also able to capture the excitonic dynamics 
over a much wider range of parameters where such sim- 
pler treatments fail. The accuracy of our variational mas- 
ter equation is verified by its excellent agreement with 
numerically converged QUAPI calculations. 

Variational transformation - As is commonjir— &2r— 
we model our QD as a two-level system with ground-state 
|0) and excited (single-exciton) state \X), separated by 
an energy of hujQ . The dot is driven by a laser of frequency 
cj;, with time-independent Rabi frequency f2, and is cou- 
pled to an acoustic phonon environment represented by 
a harmonic oscillator bath of frequencies wt, coupling 
strengths <?k, and creation operators 6^. In a frame ro- 
tating at frequency ujq, and after a rotating- wave approx- 
imation on the driving term, the total Hamiltonian reads 
(setting h = 1) 



\X)(X\ + % X +J2^bl^ (1) 



where S — ujq — uji is the detuning of the laser from the 
excitonic transition energy, (Q,/2)cr x = (n/2)(|0)(X| + 
|A}(0|) drives coherent excitonic oscillations, and an ir- 
relevant term proportional to the identity has been ne- 
glected. 
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To achieve our aim of constructing a master equation 
that is valid over as wide a range of parameters as possi- 
ble, we employ a technique similar to that originally de- 
veloped by Silbey and Harris to study the ground state of 
the spin-boson models We thus proceed by applying a 
unitary transformation to H that displaces the bath oscil- 
lators according to the state of the QD excitonic system, 
by an amount that is determined by a set {/k} of varia- 
tional parameters which we shall find through free energy 
minimization arguments. Given the particular form of 
Eq. O, we define Hy = e v He- y , where e ±v = |0)(0| + 
l*><*lllk^(±«k), with D(±Ok) = exp[±a k (^ - b k )} 
a displacement operator and = /k/wk assumed to be 
real. Notice that if /k is fixed equal to gu (for all modes) 
the transformation is equivalent to the full polaron dis- 
placement but here we explicitly give ourselves 
more freedom in choosing each /k such that we obtain 
more accurate dynamics. The idea is that by moving into 
a representation in which the displacement of the bath 
in response to the system is already accounted for, and 
by allowing some freedom to optimise the transforma- 
tion through the set {/ k }, it may be possible to identify 
a small interaction term in Hy (to be treated as a per- 
turbation) even if there is no obvious small term in the 
original Hamiltonian H. 

After transformation, the Hamiltonian reads 



R e 
Hy = —t + -a z +H E 



a.B_) + \X)(X\B z , 
(2) 



where H B = Y,k^Kb*> and a z = \X)(X\ - 
The detuning is now shifted by an amount dependent 
upon the variational parameters: e = S + R, where 
R = Ek w kVk(/k - 2p k ), which in the limit / k ->■ y k 
(full polaron displacement) becomes the polaron shift, 
R — J^k^k 5k- Note that since the term pro- 
portional to the identity in Eq. ([2|) also depends on 
the variational parameters, it cannot be neglected in 
the free energy minimisation step. The bath opera- 
tors are now given by B z = Xk(#k - /k)(& k + &k) and 
B± = ]^[ k Z3(±ak). It is straightforward to show that 
(B z ) Hb = Tt{B a exp[-pH B ]}/Tr{exp[-PH B ]} = 0, for 
inverse temperature (3 = l/kgT, but the same is not 
necessarily true for B±. We thus define B = (B±)h b , 
subtract this factor from the relevant interaction terms 
in Eq. ((2|), and add them instead to the system Hamil- 
tonian. We may then split the transformed Hamiltonian 
as Hy = Hq + Hj with the free Hamiltonian chosen to 
be H = Hs + Hb , where 



Re f2 r 
H s = —1 + -a z + —a x , 



(3) 



with Q r = Bil, while the interaction Hamiltonian be- 
comes 

Hi = Q(b x * x + B v * v )+\X)(X\B Z) (4) 

written in terms of the Hermitian combinations B x = 
(1/2) (.8+ + B_ - 2B) and B y = (l/2i)(B_ - B+). 



Importantly, in splitting Hy in the way outlined above, 
we have removed the term (Q r /2)a x from Hi to in- 
clude it instead in Hs- This has accomplished two 
things: (i) it ensures that (Hi)h = 0, which sim- 
plifies the form of the master equation to be derived 
below; (ii) more importantly, the term now appears 
as a bath-renormalised driving in Eq. ([3]) and is there- 
fore treated non-perturbatively in the subsequent for- 
malism. This ensures that our theory can capture co- 
herent exciton dynamics at the correct phonon-dressed 
energy scale f2 r . For a stationary bath state pb = 
exp[— (3H B )/Tr{exp[— /3Hb]} we find the renormalisation 
factor B = exp[-(l/2) £ k a k coth(/3w k /2)]. 

Free energy minimisation - We find the set of varia- 
tional parameters {/k} by imposing that the free energy 
associated with the transformed Hamiltonian be min- 
imised.— In doing so, we are attempting to find the best 
approximate diagonalization of H given the restricted 
form of the transformation e v He~ v , and thus expect 
contributions from the interaction Hamiltonian Hi to be 
minimised as well. In effect, the variational parameters 
allow us to tune which parts of the total Hamiltonian are 
treated exactly, and which are treated in a perturbative 
manner. 

To proceed, we compute the Feynman-Bogoliubov up- 
per bound on the free energy, given by-i^ 



A B = -iln(Tr{e- 



-PH 



}) + {Hi) Hq +O{{H 2 i) H0 ) 1 (5) 



and related to the true free energy, A, via the inequality 
Ab > A. By construction, the second term in Eq. JSJ) is 
equal to zero. Neglecting higher order terms then leads 
to the minimization condition, 



9k (l-ftanh(/V2)) 



£ tanh(/V2)(l - ^ coth(/3 Wk /2) 



1 



(6) 



which self-consistently determines the set of variational 
parameters {/k}- Here, r\ = \J t 2 + fl^ is the system 
eigenstate energy splitting in the variational frame. It is 
now insightful to look at the form of / k in two opposite 
limits (assuming, for simplicity, that we can set the de- 
tuning e to be zero): (i) when f2 << Wk we find /k ~ ffk, 
and we recover the full polaron transformation. In this 
situation, the driving is weak enough that the bath oscil- 
lators can follow the excitonic motion and become fully 
displaced when the system is in its excited state, as de- 
termined by the form of coupling in H; (ii) on the other 
hand, when Q >> cjk, we find that /k becomes very 
small, so that there is barely any displacement due to the 
transformation. In this case, the excitonic oscillations are 
too fast for the relevant bath modes to follow, and their 
displacements are thus suppressed. We shall see that 
this has important physical consequences in the context 
of driven QDs as it gives rise to a reduction in phonon- 
induced damping at large driving strengths within the 
variational theory, as originally predicted using a Feyn- 
man integral approach^ 
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Writing /k = gkF(u^) and denning the spectral den- 
sity J{uj) — J2k 1 3k 1 2 — w k), we find the following 
self-consistent forms for B and R in the continuum limit: 



B = cxp 



R 



(1/2) / J(u)F(u 
Jo 



2 coth(/3w/2)da; 



J{lj)F{lj)uj- 1 (F{lj) - 2)dw. 



(7) 
(8) 



For deformation potential coupling of the exciton state 
to acoustic phonons we can take a super-Ohmic spectral 
density of the form J(uj) = acJ^e^^'" ' ,ii2il£ where 
a captures the strength of the exciton-phonon coupling 
and uj c provides a natural high frequency cut-off, which 
is proportional to the inverse of the carrier localization 
length in the QD. 16 Throughout this paper we use the 
experimentally determined values of a = 0.027 ps 2 and 
lu c = 2.2 ps -1 ; 2 In general, for given values of a, uj c , ft, 
S and /?, Eqs. ([7]) and (JSJ must be solved simultaneously, 
which is straightforward numerically. 

Master equation - To describe the dynamics of the 
driven QD under the influence of the phonon environ- 
ment, we employ the time-convolutionless projection op- 
erator technique^ to obtain a rigorous master equation 
within the variationally-transformed representation. We 
consider the QD to be initialized in its ground state, 
for which the appropriate state of the bath is a thermal 
equilibrium state (since there is no system-bath coupling 
through the QD ground-state). We therefore write the 
initial total density operator as x(0) = |0)(0| ® p B , from 
which we find e v x(0)e~ V = x(0); our initial state is un- 
affected by the transformation. By choosing p B as our 
reference state and truncating at second-order in Hi, we 
thus find a homogeneous interaction picture master equa- 
tion 



dpv 
dt 



dsTr B [H I (t),[H I {s),p v (t)p B } 



(9) 



where pv{t) = Tr B (e v x{t)e~ v ) describes the QD exci- 
tonic degrees of freedom in the variational frame, with 
Tr b denoting a trace over the environment (phonon) de- 
grees of freedom, andO(i) = e iH t Oe ~iH t xj s i n g E q . 
and returning to the Schrodinger picture, we then obtain 



dpv 
dt 



■--i[H s ,pv(t)] 

• \ YlY,^MAi,A j {io)p v {t) - p v {t)A){w)] 

^EE s a ( w < *) ' A i W (*) + pv (*) A \ HI > 

ij id 

(10) 



where {i,j} G {1,2,3} and w E {0, ±77}. We define 
A\ = <r x , A2 — o~y and A3 = (1/2) (J + <r z ), while 
A 1 (0) = erin20(|+X+| - |-)(-|), A x {n) = cos20|-)<+|, 
A 2 (0) = 0, A 2 (r,) = M0) = cos 2 0|+)(+| + 

sin 2 I — ){ — and A^{rj) = — sin#cos# |— )(+|, defined 




FIG. 1: Exciton population dynamics calculated from the 
variational master equation (grey solid curves), weak-coupling 
theory (red dashed curves), polaron theory (blue dotted 
curves), and the QUAPI (green points), for T — 50 K. 
Here, Q, = ty/6 ps -1 (top left), Q — tt/4 ps -1 (top right), 
fl = 7r/2 ps -1 (bottom left) and Q = -k ps -1 (bottom right). 
Spectral density parameters: a = 0.027 ps 2 , lj c — 2.2 ps -1 . 



in terms of the eigenstates of H$, satisfying Hs |±) = 
(l/2)(R±rj)\±). In all cases A^u) = A\{-u), The 
angle = (1/2) arctan(f2 r /e) characterises the tilt of the 
system eigenstates away from the x-axis in the variational 
frame. 

Eq. (flUf contains the rates and energy shifts 7-y (ui,i) — 
2Rc[Kij(uj,t)] and Sij(ui,t) = lm[Kij(ui,t)], respectively, 
defined in terms of the response functions 



Kij(u,t) = [ A^e^Mr, 
Jo 



(11) 



which themselves depend on the bath correlation func- 
tions Ay(r) = Tr(Bi(T)Bj(0)p B ). Here we label B 1 = 
{fl/2)B x , B 2 = {Q/2)B y , and B 3 = B z . The bath cor- 
relation functions are most easily calculated by utilising 
the coherent state representation of the bath density op- 
erator^ We find An (t) = (fJ 2 /8)(e^ r > +e"^ T ) - 2) and 
A 22 (r) = (f7 2 /8)(e* (r) -e-*( T )), with phonon propagator 



3 ° dw ^F( W ) 2 G + (r), 



(12) 



defined in terms of G±(r) = (n(uj) + l)e~ lUT ± n(oj)e iuT , 
with n(uj) = {eP^ — 1) _1 the occupation number, while 



A 33 (r)=/ d W J( W )(l-F( W )) 2 G+(r), (13) 
Jo 

A 32 (r)=% [°° du>^-F(u)(l- F(u))iG_( T ), (14) 



with A 32 (t) = -A 23 (t), and A i2 (r) = A 2i (r) = 
A 13 (t)=A 31 (t) = 0. 

Dynamics and comparison with the path integral ap- 
proach - In order to clearly demonstrate the superiority 
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FIG. 2: Exciton population as a function of pulse area 
(O = Qt) for pulse durations of r = 2.5 ps (top) and r = 10 ps 
(bottom), and for T = 10 K (left) and T = 50 K (right). Vari- 
ational master equations results are shown by solid lines with 
the corresponding QUAPI results shown with green points. 
Spectral density parameters: a — 0.027 ps 2 , co c = 2.2 ps -1 . 

of the variational master equation over either the weak- 
coupling or polaron forms, we shall now compare its pre- 
dictions with numerically accurate results calculated us- 
ing the well-established QUAPI technique^ In all plots 
shown we have checked for numerical convergence of the 
QUAPI results. 

In Fig. [T] we plot the excited state population dynamics 
(pxx) of the QD as a function of time^ generated from 
the present variational master equation, the polaron mas- 
ter equation of Ref. the weak-coupling master equa- 
tion of Ref. HI, and our QUAPI benchmarks. For the 
realistic parameters chosen here, the variational theory 
is in excellent agreement with the QUAPI results across 
the full range of driving strengths. Importantly, neither 
the polaron nor weak-coupling theories can match the 
QUAPI results across this entire range. The polaron 
theory fails, as anticipated in Ref. |5j, at strong driv- 
ing, n > uj c . Perhaps surprisingly, the weak-coupling 
theory works well when driving beyond the cut-off fre- 
quency, even though it fails at weaker driving strengths 
due to the relatively large temperature of T = 50 K. In 
the latter regime, phonon-induced driving renormalisa- 
tion and multi-phonon effects, which the weak-coupling 
theory cannot capture, are especially important. 

In order to interpret the behaviour seen as VL is varied 
in Fig. [T] it is instructive to consider the appropriate lim- 
iting forms of the master equation (|10[) . determined by 
the bath correlation functions A^-(t). For weak driving 
conditions (though not weak system-bath coupling), we 
saw earlier that the variational transformation reduces 
approximately to performing the full polaron transfor- 
mation on H, which corresponds to the limit F(lu) —> 1 
(for all modes). In this limit, only An(r) and A22(t) 
survive and we recover the polaron master equation pre- 



sented in Ref. H, where now the perturbative expansion 
corresponds to one in the phonon-dressed driving term. 
On the other hand, when driving beyond the cutoff fre- 
quency, f2 > u> c , essentially all relevant bath oscillators 
become sluggish and they are now unable to dress the 
system, so the full polaron transformation is no longer 
appropriate. As outlined before, the correct limit is now 
F(lu) — > 0, such that effectively no transformation is ap- 
plied to the Hamiltonian and our perturbative expansion 
now corresponds to one in the original exciton-phonon 
coupling strength. Setting F(lu) = we find that only 
A33 (t) survives and our master equation reduces to the 
weak-coupling form presented in Ref. HI- In general, how- 
ever, the minimisation condition given by Eq. ([6]) corre- 
sponds to neither of these limiting cases, and F(ui) acts 
to optimise the transformation for each phonon mode 
depending on the specific Hamiltonian parameters. The 
variational master equation (|10[) therefore describes an 
intricate, yet physically motivated mixture of the weak- 
coupling and polaron approximations. 

To further confirm the versatility of the variational 
master equation, in Fig. [2] we plot the common exper- 
imental scenario of excitonic population measured as a 
function of pulse area, = fir, 1,2 for various pulse- 
widths r. Again, for all cases the variational theory 
matches the QUAPI results very closely, and in particular 
correctly predicts a recovery of Rabi oscillation amplitude 
for larger pulse areas 

Summary - We have derived a variationally-optimized 
master equation for a two-level system interacting with 
a bosonic environment, and applied it to study the dy- 
namics of phonon-damped exciton Rabi rotations in a 
laser-driven QD. By comparison with the numerically ac- 
curate QUAPI approach, we have shown that our varia- 
tional master equation can give quantitatively accurate 
results in difficult parameter regimes, where simpler mas- 
ter equation treatments fail. More generally, our results 
imply that, for super-Ohmic environments at least, the 
variational master equation can give accurate results for 
the dynamics of dissipative systems in distinct param- 
eter regimes (i.e. where the bath displacement does or 
does not play a strong role) and can reliably interpolate 
between them. Our approach is attractive not only for 
its relative simplicity, but also for the physical insight 
it gives to the system under study. In the present con- 
text, we see how the reduced damping observed using the 
Feynman integral at large pulse areas 3 - naturally arises in 
the variational theory due to the sluggish nature of the 
bath oscillators in this regime. 
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